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Abstract 

00 . 

The formation of vortex loops (global cosmic strings) in an 0(2) 
(—*) linear sigma model in three spatial dimensions is analyzed numerically. 

For over-damped Langevin dynamics we find that defect production 
is suppressed by an interaction between correlated domains that re- 
duces the effective spatial variation of the phase of the order field. 
The degree of suppression is sensitive to the quench rate. A detailed 
description of the numerical methods used to analyze the model is also 
<-h ■ reported. 
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Introduction 

The mechanism by which topological defects form in continuous phase transi- 
tions has been a fascinating field of research for many years with implications 
for astroparticle cosmology, particle physics and condensed matter systems 
|l]-{5|. On the cosmological side it is very likely that the universe underwent 
a sequence of symmetry-breaking phase transitions during its expansion and 
subsequent cooling after the Big-Bang. If the disordered phase has symme- 
try G and the ordered phase a lower symmetry H then the manifold A4 of 
possible vacua is given by the coset space G/H. In many cases the vac- 
uum manifold M. has some non-trivial topology, allowing for the appearance 
of topological defects corresponding to the non-vanishing homotopy classes 
of A4. In realistic (finite quench rate) continuous phase transitions, critical 
slowing down then implies that there must be some time during the transition 
at which the intrinsic ordering dynamics becomes too sluggish to keep pace 
with the quench. It is then inevitable that some neighboring domains will 
form with an orientation that produces topological defects as they coalesce. 
On the condensed matter side, where phase transitions occur in acces- 
sible and reproducible laboratory conditions, the production of topological 
defects is familiar and observed in, for example, ferromagnets 0, liquid crys- 
tals 0-gT| and superfluids (both A He (T2|-[T3 and 3 He |5],[16j). Q Most of 



the work on defects in condensed matter systems has focused on either the 
classification of defects || |17|-|19[1 , or the coarsening dynamics governing the 
late-time evolution of the defect density [£§. But it is also of considerable 
interest to determine the precise mechanism by which defects are produced 
and to determine the defect density at formation || ^0[ ^l|] . This is the prob- 
lem tackled in this paper for the case of the linear 0(2) sigma model, for 
which M. is the circle S . We examine numerically the process of defect 
production in fixed-rate quenches through the formation, interaction and co- 
alescence of interacting domains with well-defined phases and determine the 
defect density at the time of production. 

We find there is considerable phase alignment of domains between their 
formation and the subsequent production of defects. The spatial variation 
of the phase is then smoother than one would obtain from assuming that 
domains are statistically independent. In other words the wandering of the 



lr rhe most recent improved experiments of the Lancaster group on the fast adiabatic 
expansion of liquid 4 -He through the superfluid A transition |14f do not see any vortex 
lines, in contrast to their earlier results Itt2|. 



phase on the ground state manifold S 1 from domain to domain is not random 
- domains of a given phase attract other domains with similar phases. This 
clearly reduces the likelihood that domain coalescence will yield a topological 
defect. 

Although classical and quantum-mechanical continuous phase transitions 
have been modelled and systematically investigated, both analytically ||, 20 

and numerically p8|-[33|], a complete understanding is still lacking. 



The remainder of the paper consists of five sections. In section 1 we 
review mechanisms for the production of topological defects in continuous 
phase transitions. In section 2 we define the problem at hand and outline 
our approach. Section 3 gives the details of the numerical analysis performed 
and section 4 gives our results. Finally we conclude and discuss some possible 
implications of our work. 

1 The Kibble-Zurek mechanism 

In the Ginzburg-Landau picture continuous phase transitions proceed through 
the growth of arbitrarily small amplitude, long-wavelength fluctuations - 
the so-called spinodal decomposition of the unstable symmetric phase to the 
true ordered phase. These spinodal modes grow exponentially in time until 
cutoff by the nonlinearities associated with interactions. When a system un- 
dergoing a continuous phase transition is quenched at a finite rate from the 
disordered to the ordered phase, distinct ordered regions of space (domains) 
will generically lie at different points on the vacuum manifold. 

In [|50j Kibble proposed a clear mechanism for topological defect formation 
in a cosmological context using a simple ordering-causality argument. The 
key idea is that spontaneous symmetry breaking will occur independently 
in causally disconnected regions of space. Suppose that order parameter is 
uniform within an ordered domain of correlation volume £ 3 and randomly 
distributed on the ground-state (vacuum) manifold Ai. Furthermore Kib- 
ble assumed that the order parameter between domains was the smoothest 
interpolation possible (the geodesic rule). In this case one can, in principle, 
compute the probability that the coalescence of say three domains leads to a 
topologically non-trivial configuration. This probability is a number of order 
one so that, roughly speaking, one defect is formed per correlation volume. 

To predict the density of defects at formation in this picture one must 
have a theory of the relevant correlation length £ at formation. One simple 



proposal |20| is that £ corresponds to the correlation length at the Ginzburg 
temperature {T G ) when ordered domains are stable to thermal fluctuations. 
Given the free energy barrier AF (see Fig. 1) between the true ground state 
and the unstable high-temperature phase for T < T c , we have 



k B T G ~ e(T G )AF(T G ) 



(1) 
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Figure 1: Schematic representation of the free energy curve showing the 
disordered (T > T c ) and ordered (T < T c ) phases for a scalar field in a A0 4 
theory. The highlighted dashed section indicates typical thermal fluctuations 
for temperatures below the Ginzburg temperature T <Tq. 

Above Tq the system can locally jump back and forth between the high 
and low temperature phases. Below Tq this process is thermally suppressed. 

This argument ignores the dynamical aspects of the phase transition and 
is likely to be inaccurate if defects form at relatively high temperatures. An 
alternative non- equilibrium approach was proposed by Kibble in a later paper 
| 2l| and elaborated by Zurek ||. Consider a continuous phase transition pro- 
ceeding with a finite quench rate. The quench may be an externally imposed 



temperature or pressure quench or result from the expansion of the universe 
in the cosmological setting. In the ordered phase a given region of space 
may attain the ground state as long as the microscopic dynamics enables it 
to relax more rapidly than the quench rate. But in a continuous transition 
critical slowing down implies that the intrinsic relaxation rate becomes arbi- 
trarily slow near the critical temperature. Thus there is a characteristic time 
or temperature at which the system cannot order sufficiently rapidly. The 
correlation length at this time provides an estimate of the maximum domain 
size giving rise to topological defects. In causal language typical domain sizes 
cannot grow faster than the speed of light and therefore can never attain the 
infinite correlation lengths associated with the critical point. 
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Figure 2: Schematic representation of the relaxation time tr(£) versus the 
symmetry breaking time TsB{t) of a system being quenched through a con- 
tinuous phase transition at a rate 1/tq. 

This picture is illustrated in Fig. ^. The two time scales in the problem 
are the relaxation time tr(£) and the symmetry breaking time tsb{^) (here 
taken to be linear). The dynamics may be divided in three stages. For 
t -C —t (stage I) the relaxation time of the system is smaller than TsB(t). The 
field can dynamically relax to equilibrium while the temperature is falling. 
When t ~ —t, the symmetry breaking time and the relaxation time become 



comparable and when —t<t<t, the situation is reversed. During stage II 
T SB(t) <C r R (t), and the time needed by the system to relax to equilibrium is 
much larger than the symmetry breaking time. The system cannot relax and 
the correlation length £ cannot grow as the critical point is approached. Its 
value is frozen until, for t > t, rss(t) > Tr(£) (stage III) where £ eventually 
decreases as expected in the ordered phase. The faster the quench rate tq, 
the smaller the angle 6 indicating the slope of T SB {t). Consequently, the 
maximum value reached by £ will be smaller and smaller, creating a higher 
defect density. The slower the quench rate, the larger 9. £ is allowed to 
grow to larger values while approaching the critical point and the density of 
defects that form will be smaller (spatially correlated regions have on average 
a larger size). Assuming that the correlation length at formation is the one 
relative to t — t, (we can safely assume that £(— t) ~ £(£)), it is easy to show 
that a power law dependence between the freeze-out correlation length £* 
and the quench rate tq of the phase transition holds 

r ~ (r Q )^, (2) 

where tr(£) ~ |s|~ M ) £, ~ |s| ~ u and e = tJTQ. Extensive analytical and numer- 
ical checks have lent support to the Kibble-Zurek mechanism (see references 
in the introduction). 

2 Interacting Ordered Domains 

In this paper we will model the ordering kinetics of a non-conserved order 
parameter by Langevin dynamics with a global 0(2) <p A Ginzburg-Landau 
free energy functional. We will not treat gauge theories pi |. The order 



parameter <p is zero in the high-temperature disordered phase and non-zero 
in the low-temperature ordered phase. For a continuous transition <f) switches 
on continuously. We model the thermal quench via a linearly time varying 
mass term with slope tq. 
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i(|V0|) 2 + \/(|0|) 
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with a potential of the form 



Vm) = \m\t)\^ + \\\^ (A>0) (4) 



For m 2 (t) > 0, \<f>\ = and for m 2 (t) < 0, |0| = ±m(t)/A 1/2 . In the broken- 
symmetry phase 

4> = pe ie (5) 

where p = \cf>\ and 9 is a phase chosen from the ground-state manifold S 1 . 

Below the critical point ordering progresses via the formation of ordered 
domains within which the phase is constant. Topological defects form as 
domains coalesce, as described in the Introduction. 

When identifiable stable domains first form they are widely separated. 
Between the time they form and the time they coalesce to produce defects 
one may expect considerable interaction to take place since we are dealing 
with a non-linear field theory. Indeed angular gradient terms in the broken 
symmetry regime of the free-energy should align phases from one domain 
to the next. This will suppress wandering of the phase on the ground-state 
manifold and lower the probability of forming defects. This is the effect we 
establish and quantify in this paper. 

3 Counting Defects: A numerical approach. 

In this section we discuss how we count topological defects numerically. We 
first define an oriented (ordered) domain on the lattice as the ensemble of spa- 
tially connected sites whose phase difference satisfies the constraint \A8\ < 9 £ , 
where we have introduced a cutoff angle 9 £ . The number of domains at a 
given time is strongly dependent on the value of 9 £ . Choosing 9 £ too small 
will result in too many domains that are unstable to thermal fluctuations. 
Taking 9 £ too large, on the other hand, means we can no longer properly 
track the spatial variation of the phase. The best compromise is achieved 
by choosing the largest possible value of 9 £ that preserves the topology of 
the system. In other words, we look for the effective domains whose coa- 
lescence closely matches the distribution of defects obtained from the full 
field. To be specific we take a Z n discretization of the circle into n slices a(k) 
(k — 0, 1, ..n — 1) of width A = 2-K/n and coarse grain the angular part of 
the field as follows 

i = e a < fc > if 9 l ea(k) (6) 

where % is the lattice site and 

Q a ^ = - (jfeA + (k + 1)A) fc = 0,1,. .n-1. (7) 



We then compute the number of defects using this coarse grained field con- 
figuration and compare it to that obtained using continuous angles. We find 
that the largest A that preserves the topology is 7r/4, corresponding to a Z$ 
discretization. The error in the defect count resulting from this discretization 
is on average smaller than 5%. 

We also address the stability of domains to thermal fluctuations. For 
this purpose we introduce a minimum domain size (cut-off) A&. The system 
is kept in contact with a heat bath at constant temperature T (T <C T c ) 
throughout the simulation. A& is chosen to be the largest spatially connected 
domain generated by thermal fluctuations in the disordered phase. We now 
have all necessary tools to explore the dynamics of the effective domains. 
The following schematic illustrates our strategy: 



Evolve the system through the phase 
transition 



Discretize the phase of the field 



Build effective ordered domains out of the 
discretized field 



Compute the number of 

defects Nj generated by the 

effective domains 



X 



Randomize the phase on each 
effective domain 



Compute number of defects 

]\[Ran anc j avera gg over 3Q 

randomized configurations 



Compare N d and Nf an 



The third stage requires more details. Since we introduced a minimum 
domain size A&, very few domains will form in the early stages of ordering. 
Only a small fraction of the lattice is ordered. One may ask whether these 
domains have the expected random distribution of phases. To test this we 
simultaneously grow each domain by adding an outer shell of sites of width 
equal to one lattice site with the same field phase. We then recursively grow 
each domain's external surface until they meet in the same region of space 
and the complete lattice is filled (see Fig. |3| for an illustration of this algo- 
rithm applied to a two-dimensional configuration). The final configuration 
obtained by this construction corresponds to freely expanding domains with 




Figure 3: Two dimensional example of the domain reconstruction procedure. 
The left image represents a discretized configuration on a 50x50 lattice (pe- 
riodic boundary conditions) with a given bubble cut-off A b , and the right 
image is its layer by layer reconstructed configuration. Each color represents 
a different phase of the field. 



fixed phases. If phases are randomly distributed across domains the number 
of defects determined before and after randomizing should be statistically 
equivalent. As the quench progresses more domains are created over the lat- 
tice. We repeat the same procedure at each time step of the evolution until 
almost all the lattice is filled with well-defined domains. A domain-domain 
interaction will be traceable by comparing the number of defects produced. 
In particular the difference between the actual number of defects and the 
number for the random distribution of phases should grow with time, reach- 
ing a maximum at the time of domain coalescence. 

To resolve strings on the lattice we follow p8 and associate a vortex to 
each lattice plaquette with a non-trivial phase winding. Strings are then 
constructed by connecting these vortices, adopting a random reconnection 
algorithm for the case of multiple strings passing through the same lattice 
cell. To deal with this ambiguity the number of strings at each time step is 
obtained by averaging over 15 differently recombined string configurations 0. 
All simulations were performed on a 300MHz Pentium II for a total compu- 
tational time of roughly 800 hours. 



2 The uncertainty in this counting is smaller when the phases are discretized. 
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Figure 4: A snap-shot of a typical configuration of the system in a 3D, 
L = 30, lattice with periodic boundary conditions. Each color represents a 
different phase of the field. Dark lines represent the string defects trapped 
between the ordered domains. 

4 Numerical simulation 

We simulated a two-component classical vector field 0(r*) on a three dimen- 
sional cubic lattice of side L = 60 with periodic boundary conditions. We 
evolve the system using a leap-frog numerical implementation of the Langevin 
equation in the over-damped regime 



^30 

r~ ~dt 



v" 



dV{(f) 



+ ff(f, t) 



(8) 



where To is a dimensionful constant, V((f>) is given in Eq.(|4|) and fj(f,t) is 
Gaussian noise with temperature T that satisfies the following constraints 



(rj a {f,t))=0 



(9) 
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and 

(Va(f[, ti)^(ri, t 2 )) = 2TT 5 ab 5(fl - r$)<J(ii - t 2 ). (10) 

It is convenient to rescale to dimensionless variables, to be used from now 
on, as follows 

x = rriQX t = m t r] = rjvX/(ml) 

T' = TX/m <p' = <pV\/mQ. (11) 

where itlq is the value of m(t) prior the quench; when the system is in the 
symmetric phase. The rescaled equation, dropping the primes, then becomes 

d t( p = V 2 - a(t)<j> -<p 3 + v (12) 

where a(t) = m 2 {t)/rriQ 2 and we have chosen Tqitlq = 1. The linear quench 

is modeled by 

( 1 for t < 

a(t) = (1 - 2^) for < t < r Q (13) 

[-1 " for t>r Q 

tq~ 1 being the quench rate. This choice of a(t) enables us to drive the system 
from a disordered phase (convex potential) to an ordered phase (sombrero-like 
potential) in a finite time. For t < the system is in thermal equilibrium 
in the disordered phase. As < t < Tq, a(t) linearly decreases until it 
changes sign passing through the critical point. When t > tq, a(t) stops 
decreasing terminating the quench into the ordered phase. The limit tq — > 
corresponds to an instantaneous quench while tq — > oo describes an adiabatic 
quench. Throughout the simulation the rescaled heat bath temperature T is 
held constant. We scanned temperatures ranging from T = 0.001 to T = 0.1. 
The results reported here are for T = 0.035. The equations were numerically 
solved using Ax = 0.5 and At = 0.1. 

In the subsequent analysis we first treat instantaneous quenches and then 
compare to the results for a slow quench. 

4.1 Instantaneous Quench: tq = 

We let the system equilibrate in the disordered phase (a(t) — 1, (0) = 0) in 
the presence of thermal noise i] at T = 0.035. We sample over 1000 ther- 
malized configurations to get the minimum ordered domain cut-off A&, as 
previously described. At this temperature and this lattice size (L=60) we 
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Figure 5: Number of defects generated by the discretized field compared to 
the random-discretized field in an instantaneous quench {tq = 0). 

find that the largest thermally generated bubble in the high temperature 
phase has a volume of Vj, = A& = 79 spatially connected sites. We then 
switch a(t) = 1 —*■ ex = — 1. This breaks the 0(2) symmetry and the sys- 
tem evolves toward its equilibrium value ((f)) = 1 following the equations of 
motion defined in Eq.[L2|. We calculate the number of defects as explained 
in section || For each simulation we average the number of defects over 50 
different randomized configurations at each time step in the dynamics. The 
number of defects so obtained is then averaged over a sample of 30 differ- 
ent configurations obtained from simulations with different random initial 
conditions. 

Fig. [| illustrates the results of the analysis. At early times the number 
of defects produced by both randomized and non-randomized configurations 
agrees within statistical accuracy. This establishes that the order parameter 
is randomly varying from one spatial domain to another. At later times the 
randomized field configurations produce more and more defects with respect 
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to its counterpart. This clearly indicates that phase alignment is occurring 
as the domains grow. The gap between the two curves widens until it reaches 
a maximum at t — t* — 12 ± 1. We associate this particular time with the 
time of defect formation. This is supported by the fact that the number 
of strings generated by the effective domains reaches its maximum value at 
this same time. The number of strings subsequently decreases as domains 
coalesce and defects decay. The actual suppression factor in the number of 
strings formed is given by 7 = 2.0 ±0.2. We would like to emphasize that this 
value should be considered as a lower bound. While counting the number of 
strings we scanned our data using different minimum string length cut-offs 
A s = 4, 6, 8, 10, 12, 15, 18, 22. The data reported above reflect the analysis 
obtained with A s = 15. Smaller values of A s slightly increase the value of 7, 
but also risk over-count defects because many of the smaller strings could be 
very short-lived. Larger values of A s do not produce a detectable change in 

7- 

4.2 Slow Quench: tq = 50 

We now consider the effects of a finite-rate quench. We first equilibrate in 
the disordered phase (cx(t) = 1, (0) = 0) in the presence of thermal noise r\ 
at T = 0.035 and then initiate the linear quench according to Eq. [T3]. The 
longer time interval over which domains interact, that between the time they 
form and the time of defect production, results in a greater mismatch between 
the number of defects formed and the idealized random scenario. Of course 
this difference will disappear in the extreme adiabatic limit as no defects are 
produced at all in this equilibrium setting. We expect therefore that there is 
some finite quench rate which maximizes the suppression of defect density. 
We explored tq = 50 as a quench rate fast enough to produce a large number 
of defects but still far from the adiabatic limit. The results are given in Fig. 
[| As previously discussed, the gap between the two curves grows with time 
until it reaches a maximum value after which it falls off. It is revealing 
to compare Fig. |5| and Fig. |6| to understand the main physical difference 
between the two experiments. It is obvious that the maximum gap between 
the two curves has widened with respect to the instantaneous quench. We 
estimate the formation time to be t* = t = 30 ± 2 and 7 = 4.1 ± 0.4. The 
suppression factor 7 basically doubles while passing from tq = to tq = 50. 
There is clearly a relationship between 7 and tq. As previously discussed 
we expect that for tq — * 00, j(tq) — > and would be very interesting to 
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Figure 6: Number of defects generated by the discretized field compared to 
the random-discretized field in a slow quench (tq = 50). 

extrapolate the value of tq that maximizes j(tq) but the extremely long 
computational time required to achieve this result is prohibitive. 

5 Conclusion and Discussion 

In this paper we have analyzed the effect of a finite quench rate on the 
density of topological defects at formation. Detailed numerical simulations 
show that the phase angle of ordered domains aligns in the interval between 
domain formation and the production of defects. The effect of this alignment 
is to reduce the number of defects that form compared to the simple Kibble 
mechanism which assumes that domains are statistically independent at the 
time of coalescence. A lower bound on this relative suppression factor is 
estimated to be 7 = 2.0 ± 0.2 for an instantaneous quench and 7 = 4.0 ± 0.4 
for a phase transition with quench time tq = 50. It would be of great 
interest to systematically vary tq to determine the quench rate with optimal 
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defect suppression. As tq increases one has to explore larger correlation 
lengths which is computationally more demanding. We hope to undertake 
this challenge in the near future. 

Our results uncover a new underlying feature of defect formation mecha- 
nism that may turn out to play an important role in the precise determination 
of the initial density of defects generated in a continuous phase transition, a 
subject of great interest from both the cosmological and condensed matter 
viewpoint. Our study is most directly applicable to over-damped condensed 
matter systems but we think that an analogous effect exists in the under- 
damped regime more appropriate to a relativistic theory. An analysis in this 
direction is underway. 
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